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The  method  of  Finite  Element  Simulation,  (Balko  and 
Mendillo,  1977)  FES,  is  a powerful  numerical  technique  that 
may  he  used  to  simulate  physical  processes  in  very 
complicated  systems.  In  fact,  FES  is  sometimes  able  to 
give  useful  quantitative  results  for  a given  physical 
system,  when  conventional  analytical  and  numerical  methods 
are  difficult  to  apply. 

I 

VIhen  a physical  problem  is  solved  analytically,  the 
approach  is  to  write  a set  of  differential  equations  that 
mathematically  represents  the  physical  model  and  then  to  use 
some  well  known  analytical  method  or  educated  guess  work  to 
arrive  at  a solution.  Often  when  the  analytical  approach  is 
unsuccessful  the  same  set  of  differential  equations  can  be 
solved  numerically  using  a proven  technique  or  a synthesis 
of  several  well  known  methods.  In  the  case  of  either  an 
analytical  or  a conventional  numerical  solution,  the 
investigator  must  usually  be  very  knowledgeable  of  the 
mathematics  involved,  and  this  is  more  important  as  the 


physical  system  under  study  becomes  more  complicated  since 
Increasing  physical  complexity  greatly  Increases  the 
mathematical  complexity. 


I 


3 

being  studied,  the  familiar  equation 

J<}t>  * - Jl,  C,' 

is  the  starting  point*  Here  C|  and  C^are  the  concentrations 
of  the  reactants  and  k is  the  reaction  rate  of  the  reaction: 

C,  + C Z » Cj  -*■  Cy 

The  analytical  approach  is  to  integrate  this  equation  or  to 
write  the  complementary  equation 

C,  • 

and  to  solve  this  system  of  coupled  first  order  equations. 
In  a typical  numerical  approach  the  system  of  two  equat.iof.:' 
is  descretized  by  letting  dC— 1 * AC  and  dt-*^t.  Then 

A = ■*  ' ^5  * A C 

Though  in  the  FES  technique  it  is  not,  in  principle, 
necessary  to  write  out  the  set  of  differential  equations 
explicitly,  this  example  is  so  simple  that  FES  is,  in 
fact,  starting  with  the  differential  equations.  It  must  be 
stressed  that  this  is  a consequence  of  the  example.  Once  the 
descrete  relation 

AC  =*  - A cy  • CA  A't 

is  arrived  at  and  initial  conditions  are  deteremined,  namely 
Cj  and  at  t=0,  the  FES  technique  can  be  applied.  That 

i 8 , choose  an  aporopriate  tine  step  At,  calculate  and 

AC,,  update  the  concentrations  Cj  and  C^so  that 

C,  (t+Atj  = C,(t)  + AC, 

C±(t+At)  = Ci(t)  + ACa 


and 
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and  proceed  to  the  next  tine  step  and  calculate  in  a like 
manner.  This  is  continued  as  Ion?  as  desired. 

The  above  example  provides  a capsule  summary  of  the  basic 
components  of  the  FES  technique.  These  are: 

1.  Write  down  the  basic  equations  relating  the 
physical  parameters  of  the  system.  It  is  not  necessary 
to  explicitly  write  down  an  overall  set  of 
differential  equations  governing  the  system. 

2.  Descretize  the  basic  equations  as  necessary. 

3.  Choose  an  appropriate  time  step.  (An  important 
point  to  be  discussed  more  fully  in  Chapter  II.) 

4.  Begin  the  simulation  from  the  initial  conditions. 

5.  Determine  the  end  point  of  the  simulation,  e.g. 
some  particular  time  or  a particular  value  of  one  of 
the  parameters  being  studied. 

This  summary  suggests  that  FES  is  most  applicable  to 
initial  value  problems.  This  is  the  case.  FES  is  very 
useful  for  problems  where  the  initial  values  of  a set  of 
parameters  are  known  and  their  evolution  in  time  is  to  be 


that  are  fundamental  in  this  FES  technique. 

1.  It  is  valid  during  an  appropriately  chosen  time 
step,  to  treat  processes  that  may  be  non-linearly 
coupled  in  a linear  fashion,  i.e.,  during  the 
simulation,  to  allow  a physical  process  to  operate  as 
though  it  were  the  only  one,  then  proceed  to  the  next 
process  and  treat  that  in  a like  manner  and  to  continue 
until  all  processes  have  been  calculated  within  one 
time  step. 

2.  It  is  possible  to  determine  an  appropriate  time 
step  (which  will  fix  the  time  grid)  and  to  determine  an 
appropriate  spatial  grid. 

3.  It  is  possible  to  monitor  the  error  in  such  a 
descretized  calculation  in  some  quantitative  way. 

The  first  of  these  assumptions  is  warranted  if  the  time 
step  and  the  spatial  grid  are  of  appropriate  size.  In  fact, 
there  are  two  ways  of  handling  the  updating  of  the 
parameters.  Each  parameter  may  be  updated  after  each 


process  within  a time  step,  or  all  parameters  may  be  updated 


only  once  after  all  processes  have  operated  during  the  time 
step.  If  the  tine  step  has  been  correctly  chosen  these  two 
methods  will  give  the  same  results.  The  importance  of  the 
choice  of  the  time  step  and  the  spatial  grid  is  now  even 
more  apparent . 

The  second  assumption  is  the  most  critical.  In  order 
to  choose  an  appropriate  time  step  and  spatial  grid  the 
specific  details  of  the  physical  system  must  be  considered. 
In  general,  it  is  possible  to  do  this  by  making  common 
sense  physical  arguments  that  will  put  upper  limits  on  the 
size  of  the  time  step  as  well  as  the  spatial  grid.  The 
selection  of  a time  step  and  a spatial  grid  will  be 
illustrated  in  the  specific  case  of  the  one  dimensional 
auroral  ionosphere  model  to  be  discussed  in  the  next 
chapter. 

The  final  assumption  is  also  important  and  it  too  must 
be  justified  in  view  of  the  sp'ecific  physical  problem  being 
studied.  The  question  of  error  will  be  considered  for  the 
specific  case  of  the  one  dimensional  model  to  be  discussed 
in  the  next  chapter. 

It  is  possible  now  to  express  the  technique  of  FES  in  a 
general  mathematical  form. 
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&M] 

*(&)/>) 


= some  volume,  area,  or  linear  interval 
centered  at  location  r in  the  spatial 
grid 


= set  of  dependent  variables  being 
studied,  j=l,2,...m 


= a physical  process  relating  various 
dependent  variables  to  one  another,  to 
various  physical  parameters,  or  to  the 
independent  variables  space  and  time; 
i=  1 , 2 , • • * n 


At  = the  chosen  calculation  time  step 
which  determines  the  time  grid 


Then : 


Where 


(r-/) 


(x-i) 
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Where 


%*£>*  m change  in  QJ  at  location  due  to 

< 7 . ?\ 

process  i / J J 


Equa t ions  . I- 1 , 1-2,  and  1-3  are  used  to  calculate  at 
all  locations  of  the  spatial  grid  for  a particular  time 
step.  When  this  has  been  completed,  the  next  time  step 
proceeds  in  an  identical  manner  and  thereby  the  evolution  of 
the  in  time  is  calculated.  At  any  particular  time 
during  the  simulation  it  is  possible  to  print  out  the  h)- 
It  is  necessary  for  the  computer  code  to  store  only  one 
point  of  the  time  grid  and  all  points  of  the  spatial  grid 


for  the  values  of 


Recapitulating  the  advantages  of  FES,  first  of  all,  and 
very  important  to  the  scientist  or  engineer,  it  is  a 
technique  that  focuses  on  the  fundamental  physical  processes 
of  the  system  under  study.  The  basic  physical  relationships 
involving  the  parameters  of  int'erest  form  the  starting  point 
of  the  simulation.  There  is  no  need  to  explicitly  write  down 
a complicated  set  of  differential  equations  that  will 
describe  the  system. 


Secondly,  since  the  emphasis  is  on  the  physical 
processes  of  the  problem,  they  remain  unobscured  by 


complicated  mathematical  details  involved  in  solving  a 
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complex  set  of  differential  equations.  This  high  visibility 
of  the  physics  makes  adding  additional  physics  to  the  model 
extremely  straight  forward. 

And  finally,  there  are  many  physical  problems  in  which 
an  adequate  model  is  untractable  to  conventional  analytical 
and  numerical  techniques  (Hastings  and  Roble,  1977).  Since 
FES  is  concerned  with  the  basic  physical  processes  it  is  not 
as  susceptible  to  being  overwhelmed  by  mathematical 
difficulties  presented  by  the  problem.  If  the  individual 
physical  processes  are  known,  the  model  can,  in  principle, 
be  made  as  complex  as  desired  and  remain  tractable  to  FES. 
However,  the  computer  time  required  to  simulate  a very 
complicated  system  by  FES  may  be  great.  This  is  the  main 
draw  back  of  the  method,  but  it  is  constantly  being 
diminished  in  importance  by  the  increasing  speed  of  today's 
computers . 
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Chapter  II 

One  Dimensional  Model  of  the  Auroral  Ionoshpere 


II-l  Summary 

A one  dimensional  model  of  the  auroral  ionosphere  has 
been  devised.  The  physical  processes  that  are  included  are: 


a)  Transport:  The  transport  of  0+  and  H+  ions  due  to 
concentration  gradients  and  gravitational  drift  is 
included.  The  equation  for  the  plasma  drift  velocity 
that  was  used  (Banks  and  Holzer,  1969): 

,,  r,f-L  Jc  ZL  J-  -J~  1 

# P[c  fSi  Ti  /n(t)  <?  e H J 


Where 


D**Diffusion  coefficient 
C«Concent ra t ion  of  the  ion 
n(e)=Electron  concentration 
H«Ion  scale  height 
Te-Electron  temperature 
Ti*  Ion  temperature 

b)  Chemistry:  The  following  chemical  reactions  are 


Included : 
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N 2 + 0+ — *N0+  + 0 
02  + 0+ — ¥02+  + 0 
N0+  + e — »N  + 0 
02+  + e— *0  + 0 
H + 0+ — * 0 + H+ 

H+  + 0— *0  + + H 


c)  Production  of  0+  ions  due  to  auroral  electron 
precipitation.  In  the  two  applications  of  the  one 
dimensional  model  (winter  day  side  auroral  cusp  and 
raid-latitude  nightime  ionospheres)  there  was  no  solar 
production.  However,  if  the  model  were  to  be  used  for 
day  time  ionosphere  simulation,  solar  production  could 
be  handled  in  the  same  way  as  auroral  production. 

The  concentrations  of  the  following  species  are 
calculated  at  each  time  step:  e,  0+,  N0+ , 02+,  H+,  0,  and 

H. 

II-2  Cell  Geometry 

A one  dimensional  model  of  the  auroral  ionosphere  is 
justified  since  at  the  auroral  latitudes  where  the 
simulation  has  been  made  the  magnetic  field  lines  are 
essentially  vertical  and  parallel  up  to  a height  of  over 


2000  km.  Since  the  plasma  is  held  within  the  field  lines,  a 
first  approximation  is  Justified  where  the  cell  geometry  is 
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a stack  of  cells  for  each  of  which  the  upper  and  lower 
boundaries  are  parallel  to  the  surface  of  the  earth,  i.e.,  a 
segment  of  a circle  drawn  with  the  center  of  the  earth  as  a 
cent  er • 

Figure  II-l  shows  the  shape  and  size  of  the  boxes.  The 
center  of  the  lowest  hox  is  at  150  km  and  the  center  of  the 
top  box  is  at  2010  km.  The  reason  for  extending  the  boxes 
to  2010  km  is  to  enable  the  establishment  of  a suitable 
boundary  condition  at  the  top  cell.  This  boundary  condition 
will  be  discussed  later. 

The  cross  sectional  area  of  these  vertical  cells  is 
arbitrary  since  the  model  is  one  dimensional.  The  vertical 
dimension  is  chosen  to  be  as  large  as  possible  in  order  to 
minimize  the  number  of  cells  and  thus  the  computer  running 
time  of  the  model.  The  limit  on  the  size  of  the  cells  is 
Imposed  by  the  plasma  scale  height 
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Figure  II-l 

One  Dimensional  Model 
Cell  Geometry 


\ 


I cross  section  arbitrary 


• / 


♦ earth's  center 
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X 


Boltzman's  constant 


^ m g ( z ) , the  gravitational  acceleration 
at  height  z 


1 


Since  the  plasma  scale  height  for  H+  Is  16  tiroes  the 
scale  height  for  0+ , 0+  limits  the  size  of  the  cells.  The 
lower  16  cells  are  smaller  than  they  need  be  based  on  the 
plasma  scale  of  0+,  but  this  done  to  assure  good  spatial 
resolution  In  the  lower  F2  region  profile  just  above  and 
below  the  peak.  Twenty  50  km  cells  are  next  from  485  up  to 
1435  km.  From  1510  km  to  the  top,  the  cell  thicknesses  are 
100  km.  The  50km  and  100km  cell  heights  are  smaller  than  the 
scale  heights  of  0+  in  these  regions  but  they  were  used, 
nevertheless,  to  assure  numerical  stability  and  to  provide 
good  spatial  resolution  in  the  0+  - H+  transition  region.  In 
principle,  each  cell  could  have  been  made  a different 
thickness,  but  this  was  not  done. 

The  neutral  scale  height  was  not  directly  considered 
when  establishing  the  cell  geometry  for  two  reasons.  First, 
in  this  model  the  neutral  atmosphere  is  time  independent 
except  for  some  minor  fluctuations  in  H and  0 concentrations 
at  high  altitudes  due  to  H+  + 0 ^ * 0+  + H reactions. 
Secondly,  where  the  neutral  scale  height  is  smallest,  in  the 
lower  F region,  the  influence  of  the  neutrals  is  felt 
indirectly  through  chemistry  which  dominates  over  plasma 
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transport  as  the  most  important  process.  The  FES  treatment 
of  chemical  processes  is  independent  of  cell  size  as  will  be 
shown  later  in  the  chapter. 

The  volume  of  a cell  i is  calculated  by 

//•  - (tfe.  + ti)  ' £ * 9 ' 1 C/m 

where 


fit 

H 

1 CM 

9 


2 


= radius  of  the  earth 

= height  above  the  earth's  surface  of 
the  center  of  cell  i 
=*  arbitrary  cell  thickness 
= the  angle  (in  radians)  subtended  at 
the  earth's  center  by  the  lower  boundary 
of  the  lowest  cell 

= the  thickness  of  cell  i in  the  z 
direct  ion 


II-3  Magnetic  Field  Dip  Angle 

For  the  simulation  of  the  auroral  cusp  there  was  no 
need  to  take  into  account  the  angle  between  the  earth's 
magnetic  field  lines  and  the  vertical  cells  since  this  angle 
was  essentially  zero.  That  is  to  say,  the  diffusion  and 
gravitational  drift  of  0+  and  H+  is  parallel  to  the  cells  in 


the  auroral  cusp  region.  However,  since  it  was  decided  to 
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simulate  the  mid-latitude  ionosphere  using  initial  profiles 
measured  at  Millstone  Hill,  the  magnetic  field  dip  angle  had 
to  be  taken  into  account.  This  was  because  at  Millstone 
Hill  the  geocentric  dip  angle  is  72  deg  and  thus  diffusion 
and  gravitational  drift  would  not  take  place  along  the  same 
axis  as  the  vertically  stacked  cells.  Figure  II-2  shows  the 
geometry  of  the  situation. 

From  Figure  II-2  it  is  seen  that: 

^ 3 


Where  Vi*  is  the  gravitational  drift  velocity  in  the  z 
direction  and 

Wh  ere  is  the  diffusion  drift  velocity  in  the  z 
direction  and  Mu  is  the  diffusion  drift  velocity  along 
the  magnetic  field  line  and  I the  geocentric  dip 
angle  of  the  magnetic  field  line. 


Using  the  well  known  relation 


! c/  A/ 
A/  dJL 


Figure  1 1-2 


Magnetic  Field  Dip  Angle 
Effect  on  Plasma  Transport 


h - Direction 
/ of  Magnetic 

/ Field  Line 

/ 

/ 


1#  Geocentric  Dip  Angle 
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Where  D-  Diffusion  coefficient 

H-Humber  density  of  a species 


The  geometry  of  Figure  II-2  implies 


0 * Jl  32 


Thus 


= S/A/X 


Hence 


/ 3 <J  AA 

Ki  ^ ^a/  SiA/  x le 


Using  the  relation 


>$■  -£ 
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Where  n=Mass  of  the  ion 

k-Boltzman's  constant 
T*»Ion  temerature,  deg  K 
g“Gravi tat ional  acceleration 


Then 


]/  . e S//Ax  3 

Kj,e  £t- 
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And 

Vp.i  ^ ^>.e-  r 4,2 

vP,i  - 

Where  £ is  the  effective  plasma  velocity  in  the  z 

* * 

direction.  Conveniently,  & S/Af  .2* can  be  regarded  as  an 
effective  diffusion  coefficient  in  the  calculation  of  plasma 
transport  due  to  diffusion  and  gravitational  drift.  This 
was  done  in  the  computer  coding. 


j 

I 

[j 

li 


II-4  Neutral  Atmosphere 

The  background  neutral  atmosphere  of  the  model 

consisted  of  0,  N2,  02,  H and  He.  The  atmosphere  used  was 

O 

the  1000  K thermopause  model  from  Banks  and  Kockarts, 
(1973).  The  temperature  profile  of  this  model  neutral 
atmosphere  was  also  used.  Neutral  and  ion  temperatures  were 
set  equal  to  these  model  values.  Electron  temperature  was 
set  equal  to  twice  the  ion  temperature. 


11-5  Transport 

The  transport  of  plasma  includes  diffusion  due  to 
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concentration  gradients  and  drift  due  to  the  gravitational 
field  of  the  earth.  The  transport  of  0+  and  H+  are  handled 
in  the  same  way  so  a single,  general  explanation  may  be 
used*  Figure  II-3  shows  two  representative  cells  and  the 
relevant  geometrical  quantities. 

Let 


&,b 


The  diffusion  coefficient  of  species 


a through  species  b. 

/it  = The  overall  diffusion  coefficient  for 


spec! es  a 


Ion  concentration  in  cell  i of 


species  a 

i = Number  of  ions  in  cell  i of  species  a 
= The  calculation  time  step 
A tjfc+i  = The  cross  sectional  area  of  the 
boundary  between  cells  i and  i+1 
4+1  = The  center  to  center  distance  between 

the  cells  i and  i+1 
Zi  - The  thickness  of  cell  i 


7c  = Electron  temperature 
Tl  = Ion  temperature 

“ Electron  concentration  in  cell  i 
Vi  **  Volume  of  cell  i 


'™-l",ul  — ^—1 
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The  overall  drift  velocity,  U4, , for  an  ion  species  a 
is  given  by  (Banks  and  Holzer,  1969) 


£Ca,i  Te  t J-. 

Si  72  Ce.  22 


'I 


Where 


(Banks  and 

Kockarts,  1973) 


/L  • 


AJ2 
-n.  3(*> 


(Rishbeth  and 

Garriott,  1969) 


/>*iA 

5(p) 


Boltzman's  constant 
mass  of  ‘ion  a 

acceleration  due  to  gravity 


When  the  boundary  between  the  i and  i+1  cells  is 
considered  in  FES  this  equation  becomes 

Kin  % J 
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In  practice,  the  diffusion  coefficient  D*  . is  replaced 

/ * 

by  an  effective  diffusion  coefficient,  1)^  , given  by 


= -sw*± 


!>%i 


where  the  SwX  factor  takes  into  account  the  geocentric  dip 
angle  I and  the  term  in  brackets  weights  and  ^W^tVl  » 
the  diffusion  coefficients  within  the  cells,  according  to 
cell  size. 


The  values  used  in  the  model  for  the  diffusion  of 
0+  and  H+  were  taken  from  Banks  and  Kockarts,  (1973).  The 
values  are  shown  in  Table  II-l. 


Plasma  movement  across  the  boundary  between  cell  i and 
cell  i+1  during  a time  step  t is  described  by 

A " U,4jt-*c+l  A Cj i +\  Ciljc 

fi  r - AJA/gj  fj  r tif  Ai,i +1 

AC^i  ~ Vi  ' W'fy'+tiH  L*ti  [/% 


JCt 


4,  Itl 


9 A t>/  n r A-h  ^*A*I 

, s ~ ltA.i+c+1  Lax  At  — 
Kti 


'A,i  -*C+l  '-QjC 


(e+&)  = C«,<  (*)  + A 

£4,1+1  (ftAt)  - Cdji+M)  c+, 


-■*) 
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Table  II-l 

Diffusion  Coefficients 


cm2/sec 


D„+  - 0 - 3-71  x loU  Ti/n(0) 

dh+  - 02°  2-57  x 10  Ti/n(02) 

dH+  - N2=  2'45  x 1016  Ti/n(N2) 

D„+  - H " 7‘80  x 1016  Ti/n(H  ) 

e © 

DH+  - 0+=  6*88  X 10?  (Ti)V2/n(0+) 

DH+  - H " 8*°°  X 1017  (Ti)J5/n(H)  (14-' 

D0+  - O 3.10  x 10  (Ti)  /n(0)  (1+Tn/T^)  ** 

D0+  . 02=  7.20  X 1015  T./n{02) 

D0+  - N2=  8,50  x 10  Ti/n(N2) 

D0+  - He=  5,60  x lo16  Ti/n(He) 

D0+  - H = 4,03  x 1()16  Ti/n  (H) 

D0+  - H+=  9,39  X 10  (Ti)  ^ /n(H+) 

T^  = ion  Temperature 
Tn  = Neutral  Temperature 
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Equations  II-l  are  of  the  fora  of  equation  1-1  which 
represented  an  update  of  a general  parameter. 


In  this  manner  the  transport  of  0+  and  H+  across  all 
relevant  cell  .boundaries  Is  calculated  and  the  respective 
concentrations  are  updated  after  each  time  step. 

II-6  Chemistry 

Table  II-2  lists  the  chemical  reactions  included  In  the 
model  as  well  as  the  corresponding  reaction  rates.  The 
rates  used  are  from  Banks  and  Kockarts,  (1973).  Since  all 
the  chemical  reactions  are  handled  in  the  same  way  a single 


general  explanation  will  be  given. 


Let 


A ftC  D 

r>t  ut  t ^ ~ The  concentrations  of  four  different 
species  within  a cell 
<=  The  reaction  rate  of  a given  reaction 
At  = The  calculation  tine  step 
Then  for  the  chemical  reaction 


At  8^  Ct-  D 
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M 


Chemical  Reactions  and  Rates 


Reactions 


Reaction  Rates  (cm  /sec) 


0+  + N2  » N0+  + N k = 2 x 10 


N0+  + e » N + 0 


k = 2 x 10 


0+  + 0. 


°2+  + e 


0+  + H 


0 + H+ 


0_+  + O k = 3.46  X 10~10  /(T  )* 
^ n 


O + O k = 2.81  x 10~5  /(Te)0,85 


0+H+  k = 2. 5 x 10-11  (T  )* 

n 


0+  + H * k = 2.3  x 10-11  (T . ) ^ 


Tg  = Electron  Temperature  (Kelvin) 


T^  = ion  Temperature 


Tn  = Neutral  Temperature 


The  following  relations  hold 

AA  - -IA8A-C 
A8  = -JLA8  At 
Afc+At)  • A(t)  t AA 

S(t+A  t)  = Bit)  t A 8 


Here,  again,  is  a physical  process  changing  the 
concentrations  of  certain  species  . It  is  handled  in  the 
usual  way,  see  eq.  1-1.  It  may  be  that  the  concentrations 
of  the  products  of  the  above  reaction,  namely  C and  D,  are 
also  of  interest.  In  this  case,  by  referring  to  the  molar 
balanced  chemical  equation  (not  shown)  the  changes  in  the 
concentrations  of  the  products  C and  D can  be  related  very 
easily  to  the  concentration  changes  AA  or  AG  . This  is 
possible  since  chemical  reactions  do  not  involve  the 
transport  of  material  between  cells  that  may  be  of  different 
sizes.  Reactions  take  place  within  cells  and  therefore  the 


volume  is  constant  making  the  calculation  of  the  change  in 
the  concentration  of  a product  very  straightforward. 


II-7  Production  of  0+  ions 


The  production  of  0+  ions  due  to  the  precipitation  of 
energetic  electrons  into  the  ionosphere  was  simulated  by  the 
computer  model.  This  was  accomplished  by  using  an  0+ 
production  rate,  due  to  precipitating  electrons,  calculated 
by  Knudsen  et  al  (1977). 


■ 


A graph  of  this  production  function  is  shown  in  Figure 
11-4.  This  function  is  appropriate  for  the  winter  day  cusp 
simulation  of  the  auroral  ionosphere,  the  results  of  which 
are  presented  in  the  next  chapter. 


The  production  function  shown  in  Figure  II-4  is  given 
in  units  of  number  of  0+  ions/cn**3-sec.  In  this  form  it  is 
natural  to  handle  0+  production  in  the  model  analogously  to 
chemistry. 

Let 


At 

i>] 


■>  Production  rate  of 
(#  /cm*  * 3-sec  ) 

••  Calculation  time  step 
■ Concentration  of  0+  ions 


0+  ions 


Concentration  of  electrons 
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Then 

&[#*)  = P At 

[*%«  ■ H * 4 W 

j[e]  - a[o+] 

H.*’  H tJleJ 

The  last  two  equations  take  account  of  the  fact  that 
the  precipitating  electrons  create  ion-electron  pairs  (0+,e) 
from  neutral  0. 

When  simulating  the  winter  auroral  ionosphere  near  the 
day  side  cusp,  production  is  not  wanted  continuously.  For 
the  simulation  run  described  in  the  next  chapter  the 
production  function  was  "turned  on  and  off"  in  an  abrupt 
step  function  manner  at  the  desired  times. 


II-8  Boundary  Conditions 

In  the  section  on  cell  geometry  it  was  noted  that  the 
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bottom  cell  was  at  a center  height  of  150  kra  and  the  top  one 
at  2010  km.  The  reason  for  this  large  range  was  to 
establish  realistic  boundary  conditions.  It  would  have  been 
possible,  In  principle,  to  use  a smaller  range  and  Impose  a 
constant  plasma  flux  at  a given  height  as  a boundary 
condition.  But  It  was  decided  not  to  do  this  as  a large 
range  was  needed  anyway  since  measurements  at  the  day  side 
cusp  were  available  up  to  1400  km  and  the  problem  of 
deciding  on  a suitable  boundary  flux  could  be  avoided  as 
explained  below. 

The  boundary  condition  at  150  km  is  handled  in  an 
Implicit  way  by  allowing  the  various  physical  processes 
(diffusion,  chemistry,  gravitational  drift,  and  0+ 
production")  to  operate  normally.  The  interaction  of  these 
processes  establishes  an  equilibrium  condition  which  very 
accurately  mimics  a lower  F region  profile.  Thus,  in  the 
case  of  the  lower  boundary  condition,  there  was  no  need  to 
impose  an  ad  hoc  constraint;  'the  aeronomic  processes  alone 
generated  a realistic  profile  and  this  fact  inspires 
confidence  in  the  accuracy  of  the  model  in  the  lower  F 
region. 

There  are  two  major  differences  between  the  bottomside 
and  the  topside  that  affect  the  boundary  conditions.  First, 


on  the  bottomside 


the  neutral  atmosphere  is  very  dense  in 


F — =====  — • u 
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comparison  to  that  of  the  topside.  Thus,  even  though  the 
plasma  gradient  is  larger  on  the  bottoraside,  the  actual  rate 
of  plasma  transfer  is  small  due  to  the  relatively  small 
diffusion  coefficient  1/N  where  D is  the  diffusion 

coefficient  and  N is  the  concentration  of  the  background 
medium;  Vj  "V  D/H^  where  Vy  is  the  gravitational  drift 
velocity  and  is  the  ion  scale  height). 


The  second  reason  is  that  on  the  bottoraside  the 
dowrward  diffusing  0+  reacts  quickly  with  N2  and  02  to  form 
N0+  and  02+.  These  molecular  ions  in  turn  d i sso c i a t ive ly 
recombine  with  electrons,  thus  there  is  no  build-up  of  0+  in 
the  lower  cells  due  to  transport. 

The  case  is  different  on  the  topside.  The  neutral 
atmosphere  is  very  tenuous,  consequently,  the  diffusion 
coefficients  are  large.  Upward  diffusion  is  therefore 
large.  Since  there  is  very  little  N2  and  02  to  augment  the 
effect  of  the  "diffusive  barrier"  reaction  0+  + R ■ ♦ R+  + 
0,  there  can  be,  in  any  model  with  a finite  number  of  cells, 
a build  up  of  0+  in  the  upper  region. 

The  solution  to  the  boundary  condition  problem  was  to 
use  many  cells.  This  alone  was  not  enough,  though,  since  the 
plasma  tube  being  mimicked  extends  to  the  other  hemishpere 
and  that  is  a very  different  geometry  from  that  being  used. 
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What  was  done  was  to  add  cells  so  that  the  center  of  the  top 
cell  (100  km  thick)  reached  2010  km  and  to  not  allow  plasma 


transport 

into 

or  out  of  the 

top  cell . 

The 

effect 

o f 

this 

is  to  allow  the  chemical  reaction  H 

+ 

o 

+ 

— - H+ 

+ 0 

to 

establish 

an 

equilibrium 

between 

the 

0 + 

and 

H+ 

concentrations 

and  to  allow 

the  top 

cell 

to  act  as 

an 

inf  ini te 

sink 

« 

or  source  as 

far  as 

plasma  transport 

is 

concerned.  This  imitates  the  plasma  tube  in  a satisfactory 
way  in  this  model  since  the  region  of  interest  is  below  1500 
km  (it  is  only  at  such  heights  and  below  that  we  have  actual 
measurements  with  which  to  compare  the  results).  By 
anchoring  the  ion  concentrations  in  this  way  (the  chemical 
equilibrium  reached  will  not  differ  greatly  from  the  initial 
concentrations)  at  a height  of  2010  km  it  is  considered  that 
the  extent  to  which  this  is  an  ad  hoc  constraint  will  have 
very  little  effect  at  heights  of  1500  km  and  below. 

II-9  Calculation  Time  Step 

In  the  first  chapter  it  was  stated  that  to  illustrate 
the  selection  of  an  appropriate  time  step  the  specific 
physical  system  must  be  taken  into  account.  In  this  section 
the  physics  incorporated  in  the  model  will  be  used  to  make 
some  common  sense  arguments  to  establish  an  upper  limit  on 
the  time  step. 
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First  consider  transport. 


Le  t 


“ concentration  in  cell  i 
£\et»  = distance  between  centers  of  cells  i 


and  i+1 

%:  *=  electron  temperature  in  cell  i 


£•  *=  ion  temperature  in  cell  i 

effective  diffusion  coefficient 


between  cells  i and  i+1 
R/  = scale  height  of  the  ion  species  at 


cell  i 

Au+\  = cross  sectional  area  of  boundary 
between  cells  i and  i+1 
» volume  of  cell  i 


Consider  that  in  cell  i at  the  beginning  of  a 
calculation  time  step  there  is  ‘an  initial  concentration  of  a 
given  ion  species.  Cell  i has  boundary  interfaces  with 
cells  i+1  and  i-1  (see  Figure  II-3).  Across  those  two 
boundaries  there  can  be  plasma  transport  as  was  described 
above.  At  the  end  of  the  time  step,  after  the 
concentrations  have  been  updated,  physical  reality  requires 


that 


CcltfAt)  5 O 


■ - — — - 


OHS 
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The  plasma  velocities  across  the  two  cell  boundaries 
are  given  by  (Banks  and  Holzer,  1969) 


,,  /)  f -L-  Jj 

~ C*iHLCc  &{,£+ 1 Ct>£  % J 


fj ..  r J).  . P 1-  ~CL' t JL  £sdl£si£±  -L 


The  amount  of  plasma  transport  In  time  is 

A^itt  #<>/  6 ^ 


Ai,i-i  ~ ~ A^i-i  Cc., 


ACi  = 


Then 


C (t+A  *)  = Cfe)  * A Cc  20  3) 


By  substituting  in  for  and  working  through  the 

algebra,  eq . IT-2  can  be  shown  to  be  equivalent  to 


^ Ai-t  t?  - tU,  .v.  £ 1 (-2T  - i) 
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Following  the  terminology  introduced  by  Balko  and 
Mendillo  (1977),  it  is  convenient  to  define  a characteristic 
parameter  called  "XMOD"  for  each  physical  process  in  the 
simulation.  Each  XMOD  relates  the  time  step  to  physical 
and/or  chemical  parameters  that  are  involved  in  the  process 
it  represents. 

The  XMOD  in  eq.  II-3  is  the  dimensionless  parameter 
that  typifies  the  transpoprt  calculation.  Since  U depends 
on  the  prevailing  physical  conditions,  and  A and  V are 
fixed  by  the  cell  geometry,  XMOD,  or  t , is  fixed  by  the 
initial  conditions  but  may  be  subject  to  change  as  the 
simulation  progresses  (see  below). 

Thus  eq.  11-3,  which  is  a consequence  of 
straightforward  physical  reasoning,  sets  an  upper  limit  on 
the  only  "elastic"  parameter  appearing  in  it,  namely  /it,  the 
calculation  time  step.  This  upper  limit  can  be  used  as  a 
starting  point  in  deterraini'ng  a suitable  ^ t for  a 
calculation.  Before  discussing  how  this  is  done,  the  XMOD 
for  chemistry  will  be  derived. 

Let 


Concentration  of  a species  in  cell  i 
A set  of  n species  with  which  C 
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reacts  In  cell  i 


Jt 


= Calculation  time  step 


Assume  that  at  time  t: 


Q(t)  > ° 


Physical  reality  demands  that 


Ciit+At)  2 O 


Furthermore : 


AQ  - - Q(t)  At  <£  ^ 


i*/ 


Therefore : 


Cjt+tt)  = Cc(t)  / ^Q 


^ >o  (zr-y) 


Since  by  assumption: 


Ci  (*)  > o 


Ill* 


Equation  II-4  becomes: 


I-  At  2 O 

fz! 


* 1 (.*-*> 


>9J  ^ 

"X/VOD 


At  < 


n\ 

£4  B> 


Equation  II-5  defines  XMOD  for  a chemical  species  in 


cell  i which  reacts  with  a set  of  chemical  species  r with 


the  reaction  rates 


m 


It  is  not  necessary  to  do  an  XMOD  analysis  for  the 


other  physical  process  in  the  model,  namely  production  of  0+ 


due  to  electron  precipitation,  since  that  process  can  only 


add  plasma  and  therefore  can  never  create  a physically 


unrealistic  condition  such  as  a negative  concentration. 


The  At  must  satisfy  the  XMOD  criteria  for  transport 


and  chemistry.  Naturally  the  smaller  of  the  two  &.'b  must 


be  chosen.  In  this  way  a first  estimate  of  an  appropriate 


At  is  made.  In  practice  it  has  been  found  that  an  XMOD  as 
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large  as  0.1  yields  a time  step  that  gives  good  results. 
This  has  been  checked  by  taking  a time  step  considerably 
smaller  and  running  a calculation.  It  has  been  found  that 
the  much  finer  time  step  does  not  produce  results 
significantly  different  from  the  "coarser"  time  step  based 
on  all  XMODs  being  0.1  or  less. 

It  must  be  emphasized  that  XMOD  may  be  a function  of 
time  and  altitude.  In  the  case  of  plasma  transport,  XMOD  is 
a function  of  the  diffusion  coefficient  whch  is  a function 
of  the  species  concentrations.  In  the  case  of  chemistry, 
XMO  D is  a direct  function  of  species  concentrations. 
Species  concentrations  are  a function  of  both  time  and 
altitude. 

Another  point  regarding  XMOD  must  be  emphasized.  Each 
of  the  Xmods  was  derived  independently  of  the  others;  i.e., 
each  derivation  assumed  that  only  one  process  was  acting  and 
therefore  that  only  one  process  was  affecting  the 
concentration  of  a particular  species.  This  is  obviously 
not  the  true  picture  in  the  ionosphere  or  in  the  one 
dimensional  model.  For  example,  in  the  one  dimensional 
model  0+  is  subject  to  the  processes  of  transport  and  four 
chemical  reactions  simultaneously.  Consequently,  a time  step 
that  satisfies  each  XMOD  individually  must  be  used 
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that  satisfies  each  XMOD  individually  must  be  used 
judiciously.  If  the  calculation  were  set  up  in  such  a way 
that  all  the  concentrations  were  updated  at  the  end  of  each 
time  step  after  all  physical  processes  had  operated,  then 
the  time  step  which  satisfied  the  most  stringent  individual 
XMOD  might  not  be  sufficiently  small  to  avoid  instability 
in  the  calculation.  This  is  because,  in  this  case,  an 
"aggregate"  XMOD  (an  XMOD  calculated  taking  into  account  the 
simultaneous  action  of  several  physical  processes)  would  be 
needed  to  establish  a suitable  time  step.  What  is  done 
instead  is  to  update  the  concentrations  after  each  physical 
process  operates.  This  has  two  advantages.  First,  it 

j 

assures  the  usefulness  and  validity  of  the  individual  XMODs 
as  indicators  of  calculation  stability  and  second,  it 
enables  the  use  of  different  time  steps  for  different 
processes.  This  ootentially  saves  computer  time  if  several 
processes  have  substantially  different  time  steps  based  on 
individual  XMODs.  Recall  that  in  the  former  scheme  all 
processes  would  have  to  be  calculated  at  the  smallest  time 
step  required  by  any  one  of  the  individual  processes.  More 
will  be  said  later  on  the  use  of  multiple  time  steps. 

When  XMOD  is  used  as  a guide  for  checking  the 
appropriateness  of  the  length  of  the  time  step  it  must  be 
calculated  for  all  processes  in  all  cells.  This  may  be  done 
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XMOD  has  been  found  in  practice  to  be  a reliable 
indicator  of  when  a calculation  is  becomming  critically 
near  overshoot;  usually  manifested  by  creating  a negative 
concentration.  This  experience  has  been  gained  not  only 
with  the  one  dimensional  model  described  here,  but  also  with 
a two  dimensional  model  used  to  simulate  the  making  of 
"ionospheric  holes"  (Balko  and  Mendillo,  1977). 

11-10  Multiple  Time  Steps 

It  should  be  pointed  out  that  a simulation  need  not  use 
only  one  time  step.  In  fact  the  one  dimensional  model 
described  here  uses  five.  The  five  time  steps  are  for 
chemistry,  0+  transpoprt  in  the  lower  zone,  0+  transport  in 
the  upper  zone,  H+  transport  in  the  lower  zone,  and  H+ 
transport  in  the  upper  zone.  The  reference  to  upper  and 
lower  zones  means  that  a different  time  step  (i.e.,  a longer 
time  step)  was  used  at  lower  altitudes  for  the  transport  of 
H+  and  0+  than  at  the  higher  altitudes,  since  at  the  lower 
altitudes  the  processes  of  diffusion  and  gravitational  drift 
are  much  slower  owing  to  the  dense  neutral  atmosphere.  The 
same  time  step  was  used  for  chemistry  at  all  heights  since 
the  computer  time  required  to  calculate  chemistry  was  quite 
small  in  comparison  to  transport,  and  the  amount  of  real 


savings  was  not  worth  the  added  complexity 
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The  reason  for  using  multiple  time  steps,  of  course,  is 
to  save  computer  time.  There  is  no  point  in  an  FES 
calculation  in  using  a very  "fine"  time  step  for  any  process 
at  any  spatial  location  when  the  overall  accuracy  of  the 
results  is  determined  by  the  process  with  the  "coarsest" 
time  step  ( by  that  is  meant  a time  step  that  just  barely 
meets  the  criterion  imposed  on  it  by  the  process's  XMOD). 
The  only  complication  introduced  by  using  multiple  time 
steps  is  that  extra  computer  coding  must  be  written  to  make 
sure  that  all  processes  at  all  spatial  locations  evolve 
together  in  synchronism. 

The  ultimate  check  on  whether  a set  of  time  steps  is 
suitable  is  to  run  the  program  with  time  steps  that  are  much 
shorter  than  required  based  on  XMODs,  and  then  to  run  it 
with  as  large  tine  steps  as  possible  ( still  consistent  with 
XMODs).  This  comparison  run  need  only  be  for  a short  time, 
certainly  not  the  entire  simulation.  If  it  is  found  that 
the  two  runs  are  in  good  agreement  then  the  longer  time 
steps  are  justified.  This  approach  is  based  on  the  assumed 
accuracy  of  the  run  with  the  very  short  time  steps.  Given 
the  algebraic  form  of  the  XMODs  (which  are  a good  indication 
of  the  sensitivity  of  the  calculation  to  the  size  of  the 
time  step)  it  is  clear  that  in  every  case  decreasing  the 
time  step  can  have  only  a beneficial  effect.  This  procedure 
merely  determines  the  point  at  which  the  law  of  diminishing 
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returns  sets  in  and  enables  the  selection  of  the  largest 
feasible  time  steps  and  thus  minimizes  computer  tine* 

In  the  simulations  done  in  this  study  the  time  steps 
used  were  of  the  order  of  one  to  two  seconds. 

11-11  Accuracy 

A1 1 numerical  schemes  including  FES  are,  of  course, 
subject  to  the  errors  that  may  be  present  in  the  physics 
that  is  put  into  the  model  (for  example  diffusion 
coefficients,  reaction  rates  and  formulae  that  may  be  only 


approximately  correct).  And  all  numerical  techniques  are 
subject  to  the  errors  introduced  by  the  finite  word  length 
of  the  computer  used.  FES  has  the  advantage  , though,  of 
not  being  subject  to  another  source  of  error  that  besets 
many  numerical  techniques  and  that  is  the  errors  inherent  in 
the  algorithms  used.  Recall  that  these  techniques  are 
explicitly  concerned  with  solving  a system  of  differential 
equations.  In  order  to  achieve  that  end  approximations  are 
often  made  in  designing  the  technique  in  the  first  place. 
FES,  being  a straightforward  simulation,  is  immune  from  such 
errors.  There  is  no  mathematical  algorithm  introducing 
error  that  must  be  taken  into  account.  The  only  numerical 
source  of  error  in  the  FES  technique  is  due  to  the  finite 
computer  word  size. 
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However,  Just  as  with  other  numerical  techniques,  the 
estimation  of  these  errors  is  very  difficult  to  make.  The 
best  way  to  minimize  error  with  FES  is  to  keep  a careful 
watch  on  the  XMODs  and  to  directly  check  the  suitability  of 
the  time  steps. by  the  comparison  method. 

Another  procedure  that  proved  useful  in  assuring  the 
accuracy  and  stability  of  the  calculation  was  to  priodically 
check,  as  the  simulation  progressed,  to  see  that  in  one  time 
iteration  the  percent  change  of  each  of  the  species 
concentrations  was  not  more  than  a couple  of  percent.  This 
proved  to  be  a very  good  test.  If  the  percent  change  for 
any  important  species  became  too  large,  a judgement  was  made 
as  to  whether  to  decrease  t. 

11-12  Computer  Code 

The  one  dimensional  model  was  written  in  Fortran  for 
use  on  Boston  University's  IBM  370  computer.  It  was  divided 
into  a main  program  and  several  subroutines.  A substantial 
amount  of  input  data  is  needed  to  establish  the  cell 
geometry,  the  neutral  atmosphere,  the  initial  ionosphere  and 
the  0+  production  function.  Figure  II-5  shows  a functional 
block  diagram  of  the  computer  code. 


When  the  program  is  run  the  main  program  first  reads  in 
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some  control  data  and  the  cell  geometry.  From  then  on  It 
functions  merely  to  call  various  subroutines  which  perform 
the  actual  calculations.  The  first  subroutine  ZATMS  is 
called  only  once.  ZATMS  reads  in  the  background  neutral 
atmosphere  a nd . t emp e r a t ur es  and  also  the  initial  ionosphere. 
Some  interpolation  is  necessary  to  determine  the 
concentrations  of  the  species  In  the  various  cells  and  this 
is  accomplished  by  a subroutine  called  INTERP  which  is  not 
shown,  however,  in  Figure  11-5. 

After  ZATMS,  the  main  program  calls  ZREACT  and  ZTRANS 
in  sequence.  ZREACT  calculates  the  chemical  reactions  and 
0+  production  due  to  electron  precipitation.  ZTRANS  takes 
care  of  plasma  transport  due  to  diffusion  and  gravitational 
drift.  Since  the  diffusion  coefficients  and  gravitational 
drift  velocities  depend  ultimately  on  the  concentrations  of 
the  various  species  these  too  must  be  calculated  each  time 
step  in  ZTRANS. 

When  ZREACT  and  ZTRANS  have  completed  one  time  step  the 
main  program  checks  to  see  if  the  end  of  the  simulation  has 
been  reached  and  also  whether  a print  out  of  the  cell 
concentrations  is  required.  If  so  subroutine  AUROUT  prints 
out  the  concentrations  of  e,  0+,  and  H+  in  each  cell.  In 
addition,  the  XMODs  for  all  the  processes:  0+  diffusion, 
H+  diffusion,  the  chemical  reactions  N2  + 0+,  02  + 0+,  N0+  + 
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e,  02+  + e,  0+  + H,  and  0 + H+  are  calculated  and  printed 
out.  Also  printed  out  are  the  percent  changes  in  the 
species  concentrations  during  the  last  time  step  before  the 
printout.  Recall  that  it  Is  important  to  monitor  the  XMODs 
and  percent  changes  to  assure  that  the  calculation  does  not 
become  unstable. 

When  the  end  of  the  simulation  is  reached  a final  print 
out  of  the  parameters  is  made  by  AUROUT  and  the  program 
stops  . 
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Chapter  III 

Dayslde  Auroral  Cusp  Simulation 

III-l  Initial  Profile 

The  initial  electron  profile  used  In  the  winter  dayslde 
cusp  simulation  was  an  average  topside  vertical  electron 
density  profile  based  on  December  1971  data  from  the  Isis  2 

topside  sounder.  The  actual  profile  was  an  average  obtained 

O 

at  75.5  Corrected  Geomagnetic  Latitude  (Chacko  and  Mendillo, 
1977).  This  profile  is  shown  in  Figure  III-l. 

O 

The  75.5  CGL  profile  is  an  electron  density  profile. 
In  order  to  use  the  ionospheric  model  an  0+  and  an  H+ 
profile  were  needed.  It  would  have  been  possible  to  assume 
an  initial  H+  profile  of  zero  at  all  heights  but  that  was 
considered  unrealistic.  Instead,  the  electron  density 
profile  just  above  the  peak  was  considered  to  be  equal  to  an 
0+  profile.  A constant  0+  scale  height  was  defined  and  used 
to  determine  0+  at  all  other  topside  heights.  The  H+  profile 
was  then  calculated  by  subtracting  the  0+  profile  from  the 
electron  profile.  Thus  the  assumption  was  made  that  only 
two  ion  species  existed  at  t*=0,  0+  and  H+  (N0+  and  02  + 
profiles  were  zero  at  all  heights  at  t«0). 

III-2  0+  Production  Function 

The  purpose  of  this  simulation  was  to  mimic  the 


50 

behavior  of  the  ionosphere  beneath  the  dayside  auroral  cusp 
in  winter.  The  model  must  take  into  account,  then,  the 
production  of  plasma  by  low  energy  auroral  electrons 
precipitating  at  the  cusp.  To  do  this  a oroduction  function 
for  0+  ion-  electron  pairs  was  used  (Knudsen  et  al , 1977). 


Only  0+ 

product  ion  was 

included 

since  in 

the  region 

of 

interest 

, the  topside 

of  the 

electron 

prof ile , 

the 

production  of  0+  predominates  by  far  over  the  production  of 
N2+  and  02+.  This  production  function,  given  in  f0+/ca3-sec, 


is  shown  in  Figure  I I I- 2. 


III-3  Results 

Two  runs  were  made  in  this  simulation.  In  one  run  no 
0+  production  was  allowed  and  in  the  other  0+  production 
occurred  from  t=0  to  t=20  minutes.  This  time  interval  was 
chosen  to  test  for  the  sensitivity  of  the  high  latitude 
ionosphere  to  precipitating  electrons  beneath  the  cusp.  Both 
runs  were  carried  out  to  t=2  hours.  Comparison  can  be  made 
to  determine  the  effects  of  the  auroral  production. 


Figures  III-3  and  ITI-4  show,  respectively,  the 
profiles  at  t=20  minutes  with  no  production  and  with 
production.  It  is  seen  that  production  of  0+  ion-electron 
pairs  has  caused  a factor  of  7 enhancement  at  the  F region 
peak.  At  800km  the  effect  is  much  less  dramatic;  a factor 
of  1.3  enhancement.  This  is  expected  due  to  the  shape  of  the 
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about  640km.  This  is  also  evident  in  the  no  production  case 
shown  in  Figure  III-5.  The  enhancement  factor  at  this  lower 
H+  peak  is  3.  The  H+  profile,  as  mentioned  above,  is 
closely  coupled  to  the  0+  profile  by  chemistry  and  this 
makes  it  a sensitive  indicator  of  upward  movement  of  0+ 
enhancement  due  to  diffusion. 

The  0+  - H+  transition  height  now  shows  a substantial 
difference  between  the  two  cases.  Originally,  at  t=0,  the 
transition  height  was  at  1200km.  Now  at  t“l  hour  in  the  no 
production  case,  the  transition  height  has  moved  downward  to 
1140km  and  in  the  case  with  production  it  has  moved  upward 
to  1240km.  This  effect  can  be  understood  in  terms  of  the 
upward  diffusion  of  0+  dominating  over  loss  (pushing  up  the 
transition  height)  and,  in  the  case  without  production,  the 
loss  of  0+  at  the  higher  altitudes  due  primarily  to  downward 
transport  which  helps  to  maintain  the  F region  profile  shape 
during  the  night. 


Figures  III-7  and  TII-8  9how,  respectively,  the 
profiles  without  production  and  with  production  at  t * 2 
hours.  The  F region  peak  enhancement  has  further  declined 
to  3.4  times.  At  800kra  the  enhancement  has  further 
increased  to  3.4  times  as  the  upward  diffusion  of  0+ 
produced  at  lower  altitudes  continues.  The  lower  H+  peak 
persists  and  has  moved  upward  toward  700km.  The  0+  - H+ 
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transition  height  without  production  Is  down  to  1100km  and 
with  production  It  Is  at  about  1250km.  At  1400km  at  t“0  the 
electron  concentration  was  4 X 1 without  production  it  has 
declined  in  2 hours  to  3X1^  with  production  it  has 
increased  to  5X10  which  implies  a 1.7  times  enhancement. 


The  initial  profile  used  in  this  simulation  was  an 

0 

average  one  at  75.5  CGL.  It  was  found  (Chacko  and  Mendillo, 

1977)  that  the  ionization  produced  by  precipitating 

O 

particles  reached  a maximum  at  77.5  CGL.  It  is  interesting 

O 

to  compare  the  simulated  profile  at  75.5  with  production  and 

O 

that  averaged  from  measurements  by  Isis  2 at  77.5.  This 

O O 

assumes  that  the  quiescent  profiles  at  75.5  and  77.5  would 
be  very  similar  without  production. 


The  question  is,  then,  can  production  account  for  the 

O o 

difference  between  75.5  and  77.5?  Figures  1II-4  and  III-8 

O 

include  the  77.5  CGL  profile.  Figure  III-4  shows  the 
simulation  after  20  minutes  of'  production.  The  lower  part 
of  this  profile  agrees  reasonably  well  with  the  77.5 
profile.  At  higher  altitudes,  however,  the  difference  is 
considerable.  Figure  III-8  shows  the  "with  production" 
simulation  (which,  recall,  had  production  only  from  t»0  to 
t«20  minutes)  at  t**2  hours.  Enough  time  has  elapsed  so  that 

upward  diffusion  of  0+  has  improved  the  agreement  with  the 

O 

77.5  profile*  at  "hither  ' **alt  f tu<Tes , though*,  **  at*  l'o*w£r 
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altitudes,  the  agreement  has  worsened.  The  slope  of  the 

O 

simulated  profile  matches  the  77.5  profile  well  at  2 hours. 
If,  instead  of  only  one  20  minute  burst  of  precipitating 
electrons,  longer  "exposure  time"  to  cusp  precipitations 

were  used,  much  better  agreement  with  the  "steady-state" 

O 

77.5  CGL  profile  might  be  obtained.  It  is  clear  that  this 
model  could  be  used  in  such  a detailed  study. 


w 
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Chapter  IV 

Night-Time  Mid-Latitude  Ionosphere  Simulations 

1V-1  Introduction 

In  the  previous  chapter  the  usefulness  of  the  FES 
technique  was  demonstrated  by  the  one  dimensional  model  of 
the  auroral  ionosphere.  In  order  to  demonstrate  the 
flexibility  of  FES  it  was  decided  to  apply  the  one 
dimensional  model  to  the  mid-latitude  ionosphere.  This 
could  be  accomplished  with  minor  changes;  principle  among 
them  was  the  inclusion  of  the  geomagnetic  dip  angle.  Two 
different  initial  profiles  were  used,  both  from  Millstone 
Hill  (Evans  and  Holt, 1977).  They  differ  from  the  initial 
profile  used  in  the  cusp  simulation  in  that,  for  each,  the 
H+  - 0+  transition  height  was  800km  rather  than  1200km. 
Both  profiles  are  described  below. 


IV-2  December  6,  19'r2  Millstone  Hill  Profile 

Figure  IV-1  shows  the  Dec‘ember  6,  1 9 72  Millstone  Hill 
profile  (Evans  and  Holt,  1977).  It  was  selected  to  be  used 
as  a starting  point  for  a night-time  mid-latitude  simulation 
principally  because  of  its  rather  low  0+  - H+  transition 
height.  The  question  to  be  investigated  was  whether  such  a 
low  transition  height  profile  would  remain  stable  during 
simulation  with  the  one  dimensional  model.  The  0+  vertical 


velocity  measurements  for  this  profile  made  with  the 


incoherrent  scatter  radar  at  Millstone  Hill  were  also 


available  (Evans  and  Holt,  1977)  so  a comparison  could  be 
made  with  the  values  generated  by  the  model* 

Figure  iy-2  shows  the  profile  after  1 hour  of 
simulation.  The  peak  density  has  dropped  by  35  percent  and 
the  altitude  of  the  peak  has  risen  from  310km  to  360km.  At 
800km  the  electron  concentration  has  also  dropped  by  about 
35  percent.  The  0+  - H+  transition  height  was,  at  t=0, 
820km,  at  t**l  hour  it  has  moved  up  to  840kra.  The  peak  of 
the  H+  profile  was,  initially,  at  about  740km,  after  1 hour 
it  has  moved  upward  to  1090km  and  has  decreased  in  magnitude 
by  about  30  percent. 

Figure  IV-3  shows  the  profile  at  t=2  hours.  With 
respect  to  the  1 hour  profile:  the  peak  electron  density 
has  decreased  15  percent,  the  electron  density  at  800kra  has 
decreased  by  25  percent,  the  H+  profile  peak  has  decreased 
20  percent  and  has  remained  at  the  same  height  (1000km). 

Over  a period  of  2 hours,  then,  the  0+  - H+  transition 
height  has  remained  essentially  constant.  This  shows  that 
the  model  is  able  to  simulate  low  transition  height 
night-time  profiles. 


The  0+  vertical  -velocity  measured  at  Mil^eone  -Hil*l* 
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(Evans  and  Holt,  1977)  at  300km,  525km,  and  600km  were 

-2000cm/sec  (negative  sign  means  toward  the  earth), 
■ -A800cm/sec  and  -7500cm/sec  respectively.  The  values 
calculated  during  the  simulation  were,  at  the  same  heights 
respectively,  -2300cm/sec,  -1900cn/scc  and  -2100cm/sec.  The 
agreement  between  the  model  and  the  measurements  is  better 
than  is  first  apparent  when  the  error  bars  for  the 
observations  are  taken  into  account:  ilOOOcm/sec  at  300km, 
1f3500cra/sec  at  525kra  and  t«T4000cm/sec  at  600km. 


Another  quantity  that  can  be  compared  with  measurement 
is  the  H+  flux  at  1000km.  Evans  and  Holt  (1977)  have 
inferred  from  measurements  at  Millstone  Hill  that  in  winter 
flux  remains  upward  after  sunset,  but  usually 
becomes  downward  after  midnight  until  sunrise...".  The 
results  from  the  simulation  with  the  one  dimensional  model 
show  that  the  H+  flux  below  the  H+  peak  (1000km  at  t=l  hour) 
is  downward  and  of  the  order  of  5X10  /cn2-sec.  Above  the 
peak  the  flux  is  upward  and  about  1.0X10^ /cra2-sec.  The  fact 
that  the  simulation  does  not  indicate  a downward  flux  at 
1000km  (the  t=0  profile  is  an  after  midnight  profile)  may  be 
due  to  the  fact  that  the  upper  boundary  condition  at  2000km 
is  simply  an  infinite  sink  or  source  of  H+  and  therefore 


does  not  adequately  represent  the  coupling  between  the  night 
and  dayside  hemispheres.  This  coupling  along  plasma  tubes 
is  responsible  for  the  maintenance  of  the  downward  H+  flux 
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at  1000km  on  winter  nights. 

IV- 3 Aug.  14,  1973  Millstone  Hill  Profile 

Figure  IV-4  shows  the  Aug.  14,  1973  Millstone  Hill 

profile  (Evans  and  Holt,  1977)  at  2200  EST.  The  profile  was 
selected  for  the  simulation  because  measured  profiles  at 
five  subsequent  times  during  the  night  were  available  (Evans 
and  Holt,  1977).  It  is  characterized  by  a rather  low  0+  - 
H+  transition  height:  about  800km. 


Figure  IV-5  shows  the  profile  at  t=l  hour.  The 
profiles  shown  by  alternating  long  and  short  dashes  were 
measured  at  Millstone  Hill.  The  height  of  the  H+  peak  is  at 
650kra  for  both  the  measurements  and  the  simulation.  The 
heights  of  the  F-region  peaks  are  also  in  good  agreement, 
the  simulation  peak  being  about  40km  higher.  This  is 
understandable  in  view  of  the  characteristic  rise  of  the 
F-region  peak  as  the  F-region  decays  during  the  night.  The 
simulation,  however,  appears  to  be  "leading",  in  time,  the 
measured  evolution  of  the  profile.  The  slower  decay  of  the 
measured  profile  is  possibly  due  to  neutral  wind  transport 
of  plasma  helping  to  maintain  the  night-time  ionosphere. 
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agreement  between  the  shapes  of  the  measured  and  simulated 
profiles  - in  particular  the  slope  of  the  Ne  profile  in  the 
top-side  F-region.  There  has  been  a dramatic  change  in  the 
0+  - H+  transition  height:  it  moved  from  800km  at  t=0  to 
almost  1300km  at  t=l  hour.  This  will  be  considered  while 
examining  Figure  IV-6  which  shows  the  profile  at  t=3  hours 
(Millstone  Hill  0100  EST  A.ug.  1 5 , 19  7 3). 

In  Figure  IV-6  the  alternately  long  and  short  dashed 
lines  represent  the  measured  Millstone  Hill  profiles.  The 
shape  agreement  remains  excellent  at  t=3  hours  as  well  as 
the  agreement  between  the  peak  heights.  The  decay  rate  of 
the  model  ionosphere  continues  to  exceed  that  of  the 
measured  - this  effect  is  most  noticeable  in  the  case  of  the 
H+  profile.  The  agreement  between  the  Ne  profile  slopes 
remains  excellent.  The  0+  - H+  transition  height  in  the 
simulation  is  about  1500km  though  the  graph  does  not  show 
this  height.  The  0+  - H+  transition  height  for  the  measured 
profile  is  at  about  900km  - an  increase  of  about  100km  since 
t=0  whereas  the  simulation  shows  an  increase  from  800km  to 
1500km  In  the  same  time  span. 

Measurements  at  Millstone  were  made  of  the  vertical  0+ 
velocity.  The  values  obtained  at  Millstone  Hill  and  in  the 
simulation  at  t»3  hours  are  shown  in  Table  IV-1.  The 


agreement  in  the  flux  figures 


is  reasonable  especially  in 
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view  of  the  fact  that  the  error  bars  in  the  Millstone  Hill 
measurements  are  large  (for  example,  at  600km  and  525km 
respectively  the  uncertainties  in  the  velocity  measurements 
were  2T2100cra/sec  and  %7  30c  m / se  c . ) . The  uncertainties  in  the 


measurements  at  600km  and  525km  are  large  enough  that  even 
the  direction  of  0+  movement  is  in  doubt. 


At  t=l  hour  the  measured  H+  vertical  flux  at  1000km  was 

7 

6X10  /cm2-sec.  At  the  same  time  and  altitude  the  simulation 

7 

produced  a value  of  2.2X10/  cm2-sec.  At  t=3  hours  the 

7 

measured  value  was  1.3X10  /cm2-sec  and  the  value  resulting 
from  the  simulation  was  also  1 . 3X 1 0^ /cm2-s ec . The  agreement 
in  the  upward  H+  flux  is  very  good,  however,  it  must  be 
noted  that  since  the  flux  is  the  product  of  the  vertical 
velocity  and  the  concentration,  the  agreement  is  somewhat 
misleading  since  there  is  a large  discrepancy  in  H+ 
concentration  at  1000km  between  the  measured  and  simulated 
values.  In  the  simulation  the  flux  is  maintained  by  a 
higher  vertical  flow  velocity  ‘and  a proportionately  smaller 
concentrat ion. 


This  discrepancy  between  the  H+  profiles  deserves 

comment.  The  measured  H+  concentration  was  nearly  constant 

X 

at  a height  of  650km  from  t=0  to  t=3  hours  (3X10  /cm3).  On 
* he  other  hand  the  H+  profile  of  the  simulation  decayed  from 
' ' ' a'  at  t -0  tc  2.4X10  /cn3  at  t **  3 hours.  For  the  H+ 
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profile  peak  to  remain  in  equilibrium  for  3 hours  requires 
that  the  processes  of  chemical  production  and  loss  of  in- 


equation 


0+  + H)  as  well  as  transport  were  in  balance.  The 


. /£  - AT  = O 


must  hold  where 

is  the  H+  chemical  production  rate 
If  / cra3-sec 

is  the  H+  loss  rate,  If  / cm-sec 
AT  is  the  change  in  H+  concentration  due 
to  transport,  f /cm3-sec 

It  is  possible  for  this  equation  to  hold  at  a given  height 
for  various  H+  concentrations.  To  accommodate  this  the 
relative  values  of  Pc  Lcand  &T  would  be  different.  This  is 
the  crux  of  the  discrepancy  concerning  the  H+  concentration. 
In  the  model,  simplifications  were  made  that  affected  the 
agreement  between  the  simulation  and  the  measured  profile: 
temperature  effects  were  ignored,  i.e.  neutral, ion,  and 
electron  temperatures  were  independent  of  time;  there  were 
no  EXB  drifts  or  neutral  wind  transport  of  plasma  etc.  In 
addition,  the  1000K  neutral  atmosphere  assumed  in  the  model 
may  have  been  too  high  for  the  solar  minimum  period  being 
considered.  A cooler  neutral  atmosphere  would 
mean , ul t ima tely , that  plasma  loss  would  proceed  more  slowly 
in  the  F2  region.  These  differences  between  the  model  and 


=L 


reality,  in  addition  to  the  uncertainties  in  the  chemical 


76 


reaction  rates  and  diffusion  coefficients,  account  for  the 
differences  between  the  measured  H+  profile  and  that  of  the 
s inula t ion . 

The  ionosphere  is  a verv  sensitive  region  of  the 
earth's  upper  atmosphere  and  thus  subtle  changes  in  the 
rates  of  physical  processes  can  make  significant 
differences.  What  has  been  suggested,  though,  is  that  the 
FES  technique  shows  promise  as  a tool  to  explore  the 
sensitivity  of  the  ionosphere  to  changes  in  the  rates  of 
these  phycsical  processes. 
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Chapter  V 


Two  Dimensional  Model 


V- 1 Introduction 


The  one  dimensional  model  used  in  the  auroral 
ionosphere  simulations  was  extended  to  two  spatial 
dimensions.  The  reason  for  this  was  twofold.  First,  it 
achieves  a more  general  model  that  enables  the  simulation  of 
phenomena  that  can  be  expected  to  exibit  latitudinal 
variation,  and  second,  there  was  a specific  application  that 
could  readily  be  considered,  namely,  to  attempt  to  simulate 
the  formation  of  the  poleward  wall  of  the  mid-latitude 
ionospheric  trough  (Mendil.lo  and  Chacko,  1977)  by  asserting 
that  it  is  caused  by  the  precipitation  of  energetic 
electrons  in  the  auroral  zone.  Chapter  VI  discusses  and 
presents  the  results  of  this  simulation. 

» 

V-2  Similarities  With  The  One  Dimensional  Model 

The  two  dimensional  model  is  a direct  extension  of  the 
one  dimensional  model  and  therefore  it  has  several 
similarities.  In  fact,  several  aspects  are  exactly  the 


same . 


78 


It  was  decided  to  not  include  the  latitudinal  variation 
of  the  neutral  atmosphere  in  the  two  dimensional  model. 
Consequently,  the  same  neutral  atmosphere  as  was  used  in  the 
one  dimensional  model  was  used  throughout  the  latitude  range 
of  the  two  dimensional. 

As  described  in  Chapter  II,  the  FES  technique  treats 
chemical  reactions  in  such  a way  that  each  cell  (as  far  as 
chemistry  is  concerned)  is  independent  of  those  around  it. 
Therefore  the  inclusion  of  additional  cells  to  achieve  two 
dimensional  coverage  does  not  at  all  change  the  model's 
handling  of  chemistry.  Furthermore,  since  the  production  of 
(0+,e)  pairs  due  to  precipitating  auroral  electrons  was 
calculated  in  the  same  way  as  the  chemical  reactions  (see 
Chapter  II)  this  process  too  is  unaffected  by  the  extension 
to  two  dimensions. 


Several  other  features  of  the  model  were  altered 
substantially  enough  to  warrant  specific  discussion  below. 


V-3  Geometry  and  Dip  Angle 

Figure  V-l  shows  the  two  d imens ional  geometry  used  in 
the  model.  Also  shown  (dotted  lines)  are  a few  of  the 
earth's  magnetic  field  lines.  It  is  seen  in  Figure  V-l  that 
the  method  of  extending  to  two  dimensions  was  to  stack 
additional  one  dimensional  geometries  side  by  side,  each 


- 1 
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with  its  own  geocentric  dip  angle  as  this  parameter  is  a 
function  of  latitude.  The  "X"  axis  of  the  cells  is  parallel 
to  a geomagnetic  meridian  and  the  "Z"  axis  has  the 
advantageous  feature  of  always  being  perpendicular  to  the 
earth's  surface.  The  thickness  of  the  cells  is  arbitrary 
and  is  shown  in  Figvire  V-l  to  be  1 centimeter.  The 
boundaries  between  cells  that  are  adjacent  in  the  "X" 
direction  are  line  segments  which,  if  extended,  would 
intersect  at  the  earth's  center.  The  upper  and  lower 
boundaries  of  each  cell  are  parallel  to  the  earth's  surface. 

At  any  given  latitude,  the  number  of  cells  in  the  "Z" 
direction  and  their  respective  sizes  are  exactly  as  in  the 

one  dimensional  model.  Nine  columns  of  cells  were  used  to 

0 0 

cover  the  range  of  latitude  from  74  to  55.  The  nine  columns 

were  centered  at  74,  6 8.5*  6 7.5*  6 6.5*  6 5.5*  64.5*  63*  61* 

O 

and  55. 

V-4  Transport 

The  transport  of  H+  and  0+  ions  due  to  gravitational 


drift  and 

concentrat ion 

gradients 

is  calculated  in 

essentially 

the 

same  way 

as  in  the  one 

dimensional  model. 

see  Chapter 

II . 

Howeve  r t 

there  are  some 

complications  that 

! 

I 


must  be  described 
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First,  since  the  model  is  two  dimensional,  latitudinal 
gradients,  i.e.  in  the  "X"  direction  shown  in  Figure  V-l, 
must  be  taken  into  account.  This  is  done  by  using  the 
formulae  presented  in  Chapter  II,  but  by  omitting  the  term 
involving  the  .ion  scale  height  (gravitational  term)  since 
there  is  no  gravity  acting  in  the  "X"  direction.  There  is, 
though,  an  additional  effect  due  to  the  fact  that  the 
magnetic  field  line,  in  general,  will  not  be  aligned  with 
either  the  "X"  or  the  " Z"  axes  of  the  two  dimensional 
geometry.  This  effect  is  schematically  represented  in 
Figure  V-2.  In  the  example  illustrated,  plasma  transport 
into  or  out  of  cell(i,j)  across  the  upper  boundary  with 
cell(i,j+l)  and  the  right  hand  boundary  with  cell(i+l,j)  is 
being  considered. 

In  the  one  dimensional  model  only  the  upper  boundary  is 
relevant  and  the  effect  of  the  dip  angle  of  the  magnetic 
field  line  is  to  reduce  the  amount  of  vertical  plasma 
transport  by  the  sii?^I)  factor  (see  Chapter  II).  In  Figure 
V-2  it  can  be  seen,  however,  that  the  vertical  drift 
velocity  due  to  gravity  and  vertical  concentration  gradients 
may  be  resolved  (after  prolection  onto  the  magnetic  field 
line  along  which  all  plasma  is  constrained  to  move  - due  to 


the  gyro  frequency  being  much  greater  than  the  collision 
frequency)  into  a vertical  and  a horizontal  component.  The 
vertical  component  is  exactly  the  same  as  calculated  in  the 
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Figure  V- 2 

Magnetic  Field  Dip  Angle 
Effect  on  Plasma  Transport 


South 

U due  to  latitudinal  concentration 

gradient 


U. 
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one  dimensional  case.  The  horizontal  component,  which  was 
ignored  in  the  one  dimensional  model,  must  be  taken  into 
account  in  the  two  dimensional.  In  Figure  V-2  the  vertical 
drift  velocity  resolves  into  a downward  and  horizontal  drift 
to  the  north  (left).  There  is  no  ambiguity  concerning  the 
movement  of  plasma  as  a result  of  the  vertical  drift 
component.  Material  moves  from  cell(i,J+l)  to  cell(i,j). 
This  is  expected  since  the  gravitational  drift  velocity  was 
from  cell(i,j+l)  to  cell(i,j)  and  (in  this  example)  the 
concentration  gradient  was  calculated  between  those  same  two 
cells.  However,  what  is  to  be  done  with  horizontal 
component  of  the  drift  velocity?  It  obviuosly  arises  from 
vertical  gravitational  drift  and  vertical  concentration 


gradients.  Its  direction,  however,  indicates  that  it  moves 
plasma  from  cell(i+l,j)  into  cell(i,j).  Interestingly, 
cell(i+l,j)  has  not,  up  to  this  point,  entered  into  the 


calculation  at  all. 


This  paradox  is  a direct  consequence  of  the  discreet 
spatial  grid.  Taken  to  the  limit  of  infinitesimally  small 
cells  one  would  find  the  drift  velocity  at  a point  in  space 
to  be  along  the  magnetic  field  line.  This  is  consistent 
with  differential  equations  governing  this  physical  system. 
As  has  been  stressed  though,  FES  i3  not  explicitly  concerned 
with  the  differential  equations.  Therefore,  a convention 


must  be  adopted  to  resolve  this  difficulty  that  is  typical 
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of  those  that  arise  when  an  essentially  continuous  process 
is  modeled  in  a discreet  manner.  It  should  he  noted,  also, 
that  any  discreet  model  would  have  a similar  difficulty  and 
that  this  current  problem  is  not,  then,  peculiar  to  FES. 

The  convention  used  is  straightforward.  Whenever  the 
transport  into  or  out  of  a given  cell(i,j)  is  considered  the 
concentration  gradients  across  only  two  boundaries  are  used; 
that  with  cell(i,j+l)  and  that  with  cell(i+l,j).  If  this  is 
done  for  all  cells(i,j)  then  it  is  easy  to  visualize  that 
all  boundaries  are  taken  into  account  and  that  when  the 
update  of  plasma  concentrations  is  made  after  having 
considered  all  cells,  all  transport  among  adjacent  cells  has 
been  included . 

This  convention  implies  that,  in  general,  (beyond  the 
specific  case  illustrated  in  Figure  V-2)  any  vertical  drift 
component  moves  plasma  across  the  boundary  between 
cells(i,j)  and  (i,j+l)  and  that  any  horizontal  drift 
component  moves  plasma  across  the  boundary  between 
cells(i,j)  and  (i+l,j).  Thus  is  preserved  the  computational 
convention  that  when  transport  into  or  out  of  cell(i,j)  is 
considered,  transport  may  occur  only  across  the  upper  and 
right  hand  boundaries.  Since  the  cell  sizes  used  in  the 


model  are  appropriate  for  adequate  spatial  resolution  this 
method  yields  good  results. 
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V-5  Boundary  Conditions 

indicated  that  one  dimensional  cell  geometries  were  stacked 
side  by  side  to  embody  the  two  dimensional  geometry.  It  was 
natural,  then,  to  use  exactly  the  same  boundary  conditions 
at  the  top  and  the  bottom  of  the  two  dimensional  model  as 
were  used  in  the  one  dimensional  (see  Chapter  II). 


The  two  dimensional  model,  however,  involves  two  more 
boundary  conditions:  those  at  the  northern  and  southern 
latitude  limits.  During  the  development  of  the  one 
dimensional  model  it  was  learned  that  care  must  be  taken  to 
assure  that,  at  the  higher  altitudes,  the  cells  at  the 
extremities  of  the  region  do  not  tend  to  "fill  up" 
appreciably  during  the  simulation.  As  already  mentioned,  in 
the  case  of  the  top  cells,  this  was  prevented  by  using  the 
same  boundary  conditions  as  in  the  one  dimensional  model. 

The  boundary  conditions  used  in  the  two  dimensional 
model  at  the  northern  and  southern  extremeties  were  the 
same.  What  was  done  was  to  not  allow  the  plasma 
concentrations  in  the  columns  of  cells  forming  the  northern 
and  southern  boundaries  to  change  in  response  to  latitudinal 
plasma  transport.  The  reason  that  this  very  simple  boundary 


condition  worked  is  that  the  latitudinal  transport  of  plasma 
at  the  mid  and  high  latitudes  (the  range  of  the  trough 
poleward  wall  simulation)  is  rather  small  due  to  the  small 
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latitudinal  concentration  gradients  and  to  the  orientation 
of  the  magnetic  field  lines.  Another  reason  is  that  the 

cells  at  the  northern  and  southern  extremities  were  made 

O 

large  (10  in  extent). 

If  the  model  were  to  be  used  to  simulate  the  ionosphere 
in  the  equatorial  region  it  might  be  necessary  to  use  a 
different  boundary  condition  at  the  latitudinal  extremes. 
This  would  be  necessitated  by  the  orientation  of  the 
magnetic  field  lines  which,  near  the  equator,  would  channel 
plasma  movements  in  latitudinal  directions. 


V-6  Calculation  Time  Step  and  XMODs 

By  extending  the  model  to  two  dimensions  an  added 
degree  of  freedom  was  given  to  the  process  of  transport.  In 
the  one  dimensional  model  transport,  of  course,  could  only 
occur  in  the  vertical  direction,  but  in  the  two  dimensional 
it  could  also  occur  horizontally.  To  answer  the  question  - 
Would  the  same  time  steos  that  were  suitable  in  the  one 
dimensional  model  also  work  in  the  two  dimensional?  - an 
XMOD  calculation  for  horizontal  transport  was  made.  As 
expected,  the  process  of  horizontal  transport  is  quite  slow 
due  to  the  smalt  latitudinal  concentration  gradients  and  to 
the  orientation  of  the  magnetic  field  lines.  Thus 


horizontal  transport  could  use  a longer  time  step  than 
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vertical  transport  if  desired.  It  was  decided,  however,  to 
use  the  sane  time  step  for  all  transport. 

The  latitudinal  variation  of  the  ionosphere  over  the 
range  considered  in  the  poleward  wall  simulation  was  snail 
enough  to  permit  calculating  XMODs  for  all  processes  along 
only  one  column  of  cells  at  a particular  latitude.  This 
saved  computer  tine  and  facilitated  the  interpretation  of 
the  XMODs  when  checking  the  stability  of  the  calculation. 
Had  there  been  a large  variation  with  latitude  of  important 
physical  parameters,  then  XMOD  calculations  would  have  been 
necessary  at  each  column  of  cells  to  make  sure  that  any 
instability  in  the  simulation  would  not  go  unnoticed . 

In  summary,  then,  the  same  time,  steps  for  all  processes 
and  at  all  spatial  locations  were  used  in  the  two 
dimensional  simulation  as  were  used  in  the  one  dimensional. 
These  were  on  the  order  of  one  to  two  seconds. 

V-7  Accuracy 

Extending  the  model  to  two  dimensions  did  not  introduce 
any  new  sources  of  error  in  the  simulation  nor  did  it  change 
any.  This  is  one  of  the  advantages  of  FES.  The  finite  word 
length  of  the  computer  is  the  primary  source  of  error 
regardless  of  the  complexity  of  the  model. 
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This  very  often  is  not  the  case  with  alternative 
methods  of  numerical  analysis  of  complex  systems.  Those 
techniques  that  are  explicitly  concerned  with  solving 
differential  equations  are  not  so  easily  extended  to 
additional  spatial  dimensions.  This  is  because,  frequently, 
entirely  different  mathematical  algorithms  must  be  developed 
or  adaoted.  When  the  algorithms  are  changed  to  accommodate 
greater  mathematical  complexity  it  is  not  uncommon  that  the 
overall  accuracy  of  the  method  will  deteriorate  beyond  that 
caused  by  the  additional  round-off  errors.  This  situation 
means  FES,  in  principle,  has  an  advantage  as  no  algorithms 
are  ch  anged  . 
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Chapter  VI 

Mid -Latitude  Trough  Poleward  Wall  Simulation 


VI-1  Introduction 

The  two  dimensional  model  was  described  in  detail  in 
the  previous  chapter.  The  simulation  of  the  mid-latitude 
trough  poleward  wall  was  selected  as  an  interesting  first 
application.  The  geometry  of  the  model  is  shown  in  Figure 
V-l  . 


VI-2  Initial  Two  Dimensional  Profile 

Figure  VI-1  shows  the  initial  profile  used  in  the 
simulation.  It  Is  based  on  the  average  ionospheric  trough 
found  by  Mendillo  and  Chacko  (1977)  during  a period  of  very 
low  magnetic  activity  in  December  of  1971.  In  order  to 
investigate  whether  particle  precipitation  in  the  auroral 
zone  could  be  responsible  for  the  shape  of  the  poleward  wall 
of  the  trough,  the  initial  profile  was  established  as  shown 

in  Figure  VI-1  with  constant*  Ne  concentrations  at  all 

O 

latitudes  northward  of  64.5  Corrected  Geomagnetic  Latitude 
(CGL) . The  simulation  employed  an  0+  ion-electron  pair 
production  rate  typical  for  the  auroral  zone  (Knudsen  et  si, 
1977)  that  is  shown  in  Figure  VI-2.  The  expectation  was 
that  the  characterisitc  shape  of  the  poleward  wall  found  by 
Mendillo  and  Chacko  (1977)  would  emerge  during  the 
simulation.  Since  the  poleward  wall  was  found  to  occur  at 
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about  64  CGL  (Mendillo  and  Chacko,1977)  which  is  about  1 - 

O 

2 equatorward  of  the  statistical  location  of  the 
equatorward  edge  of  the  diffuse  aurora,  the  production 
function  was  "allowed  to  operate"  only  in  three  columns  of 

cells  shown  in  Figure  V-l  (all  other  cells  had  no 

o 

production).  These  three  columns  spanned  3 in  latitude  from 

O • 

66  CGL  to  69  CGL;  the  production  of  0+  ion-electron  pairs 
in  this  region  mimicked  the  effect  of  the  diffuse  auroral 
pr ecipitat ion. 

Vl-3  Results 

At  t»0  the  production  of  0+  ion-electron  pairs  began. 
Thfe  production  function  used  to  mimic  the  diffuse  aurora  was 
"turned  on"  for  the  initial  30  minutes  of  the  simulation  and 
"turned  off"  for  the  remaining  30  minutes.  This  coarse  time 
modulation  of  the  production  was  selected  because  it  was 
simple  to  impose  and  adequate  to  get  an  idea  of  the  effect 
of  this  production  on  the  development  of  the  poleward  wall. 

Figure  Vl-3  shows  the  profile  at  t-30  minutes  (at  the 
time  when  production  was  turned  off).  Figure  VI-5  shows  an 
average  steady-state  poleward  wall  profile  (Mendillo  and 
Chacko,  1977)  for  a quiet  magnetic  period.  The  enhancement 
at  1^^  and  450km  in  the  simulation  is  considerably  greater 
than  that  appearing  in  Figure  VI-5.  At  550km  the  agreement 
between  the  simulation  and  Figure  VI-5  is  excellent;  above 


this  height  the  development  of  the  wall  is  less  than  the 
steady-state  profile.  This  may  he  explained  by  the  fact 
that  the  production  function  (see  Figure  VI-2)  is  relatively 
large  from  about  150kra  to  about  450kra;  outside  this  range 
the  production . of  0+  is  rather  insignificant.  Therefore  the 
formation  of  the  poleward  wall  at  the  higher  altitudes  is 
dependent  on  the  process  of  upward  diffusion  from  the  lower 
altitudes  where  production  is  substantial.  At  t=30  minutes 
the  large  amount  of  0+  created  at  the  lower  altitudes  has 
not  had  an  opportunity  to  diffuse  upward. 

Figure  VI-4  shows  the  simulation  profile  at  t=*60 
minutes  ( from  t=30  minutes  to  t=60  minutes  there  was  no 
production  of  0+) . Upward  diffusion  has  been  smoothing  out 
the  vertical  Ne  gradient.  The  poleward  wall  is  now  better 
developed  at  650km  through  950km.  At  1400km,  however,  the 
profile  has  decayed  more  than  it  has  been  enhanced  by  upward 
diffusion;  therefore,  little  of  the  poleward  wall  shows  up 
at  this  height  and  this  is  the  ‘cause  for  the  poor  comparison 
with  the  steady-state  profile  of  Figure  VI-5. 

At  1^  , 450km  and  550km  the  simulation  profile  has  a 
higher  peak  density  at  the  "top  of  the  wall"  than  does  the 
steady-state  profile  - 60,  65  and  36  percent  respectively. 
This  is  due  to  the  fact  that  the  simulation  profile  at  t“60 
minutes  is  not  in  equilibrium  and,  of  course. 
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production  function  of  Figure  VI-2  that  was  turned  on  for  30 
minutes  and  then  off  for  30  minutes  is  probably  not  the  best 
approximation  to  the  real  steady-state  condition. 
Nevertheless,  the  agreement  between  the  average  steady-state 
profile  and  the  simulation  is  rather  good  especially  with 
respect  to  the  general  slope  of  the  poleward  wall  at  the 
various  altitudes. 

Figure  VI-6  and  VI-7  show  the  bottomside  F region 
profile  at  t'  = 0 and  t = 60  minutes  respectively.  There  are  no 
measurements  with  which  to  compare  as  the  Mendillo  and 
Chacko  analysis  (1<177)  was  based  on  Isis  2 topside  sounder 
data.  The  initial  bottomside  profile  shown  in  Figure  VI-6 
was  arbitrarily  chosen. 

The  fact  that  the  poleward  wall  of  the  mid-latitude 
ionospheric  trough  is  attributable  to  the  precipitation 
associated  with  the  diffuse  aurora  is  well  known  and  this 
finite  element  simulation  supports  the  plausibility  of  such 
an  explanation.  The  hope  is  that  the  flexibility  of  FES 
(particularly  the  ease  with  which  a model  may  be  extended  to 
another  spatial  dimension)  has  been  further  demonstrated. 
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Sunrna  ry 

A two  (I  t me  ns  lon.i  l conputer  model  of  the  night-time 
ionosphere  has  been  developed  to  simulate  atmospheric 
processes  acting  at  mid  and  high  latitude  regions.  Since 
the  temporal  scales  of  the  effects  considered  were  much  less 
than  a day,  an  "initial  value  problem"  apnroach  was  adopted. 
The  computational  method  used  was  a numerical  technique 
whereby  the  finite  clement  simulation  (FES)  approach  was 
applied  to  physical  and  chemical  processes,  rather  than  an 
explicit  solving  of  the  set  of  differential  equations 
governing  the  ionospheric  system. 

The  following  physical  processes  were  included  in  the 
model:  transport  of  0+  and  ii+  ions  due  to  am  bipolar 
diffusion  caused  by  concentration  gradients  and 
gravitational  drift  (transport  of  N0+  and  02+  was  not 

treated),  the  chemical  reactions  0+  + 112 ► 'N  + N0+,  H0+  + 

e ►N  + 0,  0+  + 02  * 02+  + 0,  02+  + e ► O + 0 , H + 0+ 

11+  + 0,  production  of  0+  ions  due  to  the  precipitation  of 

auroral  electrons. 

A neutral  atmosphere  consisting  of  0,  02,  N?  , He,  and  H 
was  input  to  the  model  as  well  as  an  initial  plasma  profile 
consisting  of  H+ , 0+,  H0+  and  02+  (>T0+  and  02  + 
concentrations  are  both  zero  initially).  During  simulations 


the  temperature  profiles  of  the  neutral  atmosphere  and 
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plasma  did  not  change  (electron  temperature  was  everywhere 
equal  to  twice  the  ion  temperature  which  was  everywhere 
equal  to  the  neutral  temperature).  Furthermore,  the 
concentration  profile  of  the  neutral  atmosphere  was 
considered  essentially  time  independent. 

The  cell  geometry  was  an  array  of  "wedge"  shaped  cells 
with  northern  and  southern  boundaries  formed  by  rays  drawn 
from  the  earth's  center.  The  upper  and  lower  boundaries  of 
each  cell  were  arcs  drawn  with  a geocentric  center.  The  two 
spatial  dimensions  covered  by  this  geometry  were  altitude 
(from  IS 0km  to  2 000km)  and  latitude  (for  the  poleward  wall 

simulation  the  latitudinal  extent  of  the  array  was  about 

o 

20) . The  north-south  axis  ot  the  array  was  colinear  with  a 
geomagnetic  meridian.  Suitable  boundary  conditions  for 
physical  processes  were  used  at  the  extremities  of  the 
array.  In  altitude  thickness  the  cells  ranged  from  20km  to 
100km  at  low  and  high  altitudes  respectively. 

When  a simulation  was  made  the  initial  profiles 
(neutral  and  plasma)  were  input  and  the  time  evolution  of 
the  concentrations  of  0+,  '!+,  N0+,  and  02+  was  calculated 
for  each  cell.  This  was  carried  out  using  a suitable  time 
step  for  each  process  and  updating  the  ton  concentrations  in 
response  to  chemical  and  transport  effects. 
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A one  dimensional  (height-dependent)  version  of  the 
code  was  used  to  simulate  the  ionosphere  beneath  the 
mid-winter  dayside  auroral  cusp  (effectively  a night-time 
ionosphere).  The  one  dimensional  code  was  also  used  to 
simulate  the  behavior  of  night-time  Millstone  Hill  profiles. 
A two  dimensional  (height  and  latitude)  simulation  was 
carried  out  For  the  mid-latitude  ionospheric  trough  and 
poleward  wall  region.  In  each  case,  the  simulations  spanned 
a 2-3  hour  time  period. 

This  work  has  shown  the  flexibility  of  the  FES 
technique  and  its  suitability  for  computer  modeling  studies 
of  ionospheric  processes  and  effects.  In  the  future, 
f r>mn er ?. t ur ? effects  and  clactrcd'/n amic  drifts  could  be  added 
in  order  to  achieve  a more  realistic  model.  beyond  these 
two  processes,  the  effects  of  neutral  wind  plasma  transport 
and  additional  chemical  reactions  might  also  be  included. 
Finally,  an  extended  model  (including  the  additional 
processes  mentioned  above)  could  be  expanded  to  3 spatial 
dimensions  to  realize  the  most  general  model  possible. 
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